
clear;
%insert "Data" file location in [...]
folder = '[...]Data';
cd(folder)
addpath(folder)

para0found = []
para1found = []

para0found_n = []
para1found_n = []

hrl0all = []
trl0all = []
hrl1all = []
trl1all = []

tch0 = []
tch1 = []
tnch0 = []
tnch1 = []

tchall0 = []
tchall1 = []

rhealthier0all = []
rhealthier1all = []
rtastier0all = []
rtastier1all = []

right0all=[]
right1all=[]


varRall=[]
Nall=[]

for i =  [13:25 27:29 31:32 34:37 39:55 57:60 62:65 66:68]
    X=readtable (strcat('Efood',  string(i), '_answers_new', '.xlsx'));
health=X.('Health')
taste=X.('Taste') 
nH=nansum(sign(abs(health)+1))
nT=nansum(sign(abs(taste)+1))
varR=((nH-1)/(nH+nT-2))*(nanstd(abs(health)))^2 + ((nT-1)/(nH+nT-2))*(nanstd(abs(taste)))^2 + 2*corr(abs(health),abs(taste),'rows','complete')*nanstd(abs(health))*nanstd(abs(taste))*((nT-1)/(nH+nT-2))*((nH-1)/(nH+nT-2))
varRall=[varRall;varR*(nH+nT)]
Nall=[Nall;nH+nT]
end

    avstdR = sqrt(sum(varRall)/sum(Nall))
   
probitf0ch = []
probitf1ch = []
probitf0nch = []
probitf1nch = []

tprobitf0ch = []
tprobitf1ch = []
tprobitf0nch = []
tprobitf1nch = []

taprobitf0ch = []
taprobitf1ch = []
taprobitf0nch = []
taprobitf1nch = []

hprobitf0ch = []
hprobitf1ch = []
hprobitf0nch = []
hprobitf1nch = []


for i = [13:25 27:29 31:32 34:37 39:55 57:60 62:65 66:68]
 
    
T0=readtable (strcat('behavior',  string(i), '_output_results0', '.xlsx'));
T1=readtable (strcat('behavior',  string(i), '_output_results1', '.xlsx'));

type0 = T0.('Type')
type1 = T1.('Type')

for r=1:90
   if type0(r,1)==1|type0(r,1)==3 
    chall0(r,1)=1
   else
    chall0(r,1)=0
   end 
end

for r=1:90
   if type1(r,1)==1|type1(r,1)==3 
    chall1(r,1)=1
   else
    chall1(r,1)=0
   end    
end

ch0 = chall0
nch0 = 1-chall0
ch0(chall0==0)=NaN
nch0(chall0==1)=NaN

ch1 = chall1
nch1 = 1-chall1
ch1(chall1==0)=NaN
nch1(chall1==1)=NaN

tch0 = [tch0, ch0]
tch1 = [tch1,ch1]
tnch0 = [tnch0, nch0]
tnch1 = [tnch1, nch1]

tchall0 = [tchall0, chall0]
tchall1 = [tchall1, chall1]

sel0=T0.('Selection') 
hrl0=(T0.('Stimuli_right_Health') - T0.('Stimuli_left_Health')) /avstdR
trl0=(T0.('Stimuli_right_Taste') - T0.('Stimuli_left_Taste')) /avstdR
sel1=T1.('Selection')
hrl1=(T1.('Stimuli_right_Health') - T1.('Stimuli_left_Health')) /avstdR
trl1=(T1.('Stimuli_right_Taste') - T1.('Stimuli_left_Taste')) /avstdR

for r=1:90
if sel0(r,1)=="right"
    right0(r,1)=1
elseif sel0(r,1)=="left"
    right0(r,1)=0
else
    right0(r,1)=NaN
end
end

for r=1:90
if sel1(r,1)=="right"
    right1(r,1)=1
elseif sel1(r,1)=="left"
    right1(r,1)=0
else
    right1(r,1)=NaN
end
end

rhealthier0 = sign(hrl0)
rhealthier1 = sign(hrl1)
rtastier0 = sign(trl0)
rtastier1 = sign(trl1)

hrl0all = [hrl0all,hrl0]
trl0all = [trl0all,trl0]

hrl1all = [hrl1all,hrl1]
trl1all = [trl1all,trl1]

rhealthier0all = [rhealthier0all,rhealthier0]
rhealthier1all = [rhealthier1all,rhealthier1]
rtastier0all = [rtastier0all, rtastier0]
rtastier1all = [rtastier1all, rtastier1]

right0all=[right0all,right0]
right1all=[right1all,right1]


mch0 = glmfit([hrl0 trl0 rhealthier0], right0.*ch0, 'binomial','Link','probit')';
mch1 = glmfit([hrl1 trl1 rhealthier1], right1.*ch1, 'binomial','Link','probit')';
mnch0 = glmfit([hrl0 trl0 rhealthier0], right0.*nch0, 'binomial','Link','probit')';
mnch1 = glmfit([hrl1 trl1 rhealthier1], right1.*nch1, 'binomial','Link','probit')';

probitf0ch = [probitf0ch;mch0]
probitf1ch = [probitf1ch;mch1]
probitf0nch = [probitf0nch;mnch0]
probitf1nch = [probitf1nch;mnch1]

if mch0(1,2)>0&mch0(1,3)>0
    mch0=mch0
elseif mch0(1,2)>0&mch0(1,3)<0
mch0 = glmfit([hrl0 rhealthier0], right0.*ch0, 'binomial','Link','probit')'
mch0 = [mch0(1,1), mch0(1,2), 0, mch0(1,3)] 
elseif mch0(1,2)<0&mch0(1,3)>0
    mch0 = glmfit([trl0 rhealthier0], right0.*ch0, 'binomial','Link','probit')'
    mch0 = [mch0(1,1), 0, mch0(1,2), mch0(1,3)] 
else
    mch0=[NaN,0,0,0]
end

if mch1(1,2)>0&mch1(1,3)>0
    mch1=mch1
elseif mch1(1,2)>0&mch1(1,3)<0
mch1 = glmfit([hrl1 rhealthier1], right1.*ch1, 'binomial','Link','probit')'
mch1 = [mch1(1,1), mch1(1,2), 0, mch1(1,3)] 
elseif mch1(1,2)<0&mch1(1,3)>0
    mch1 = glmfit([trl1 rhealthier1], right1.*ch1, 'binomial','Link','probit')'
    mch1 = [mch1(1,1), 0, mch1(1,2), mch1(1,3)] 
else
    mch1=[NaN,0,0,0]
end

if mnch0(1,2)>0&mnch0(1,3)>0
    mnch0=mnch0
elseif mnch0(1,2)>0&mnch0(1,3)<0
mnch0 = glmfit([hrl0 rhealthier0], right0.*nch0, 'binomial','Link','probit')'
mnch0 = [mnch0(1,1), mnch0(1,2), 0, mnch0(1,3)] 
elseif mnch0(1,2)<0&mnch0(1,3)>0
    mnch0 = glmfit([trl0 rhealthier0], right0.*nch0, 'binomial','Link','probit')'
    mnch0= [mnch0(1,1), 0, mnch0(1,2), mnch0(1,3)] 
else
    mnch0=[NaN,0,0,0]
end

if mnch1(1,2)>0&mnch1(1,3)>0
    mnch1=mnch1
elseif mnch1(1,2)>0&mnch1(1,3)<0
mnch1 = glmfit([hrl1 rhealthier1], right1.*nch1, 'binomial','Link','probit')'
mnch1 = [mnch1(1,1), mnch1(1,2), 0, mnch1(1,3)] 
elseif mnch1(1,2)<0&mnch1(1,3)>0
    mnch1 = glmfit([trl1 rhealthier1], right1.*nch1, 'binomial','Link','probit')'
    mnch1= [mnch1(1,1), 0, mnch1(1,2),  mnch1(1,3)] 
else
    mnch1=[NaN,0,0,0]
end

%final probit coefficients (with positive Hrl and Trl coefficients)
tprobitf0ch = [tprobitf0ch;mch0]
tprobitf1ch = [tprobitf1ch;mch1]
tprobitf0nch = [tprobitf0nch;mnch0]
tprobitf1nch = [tprobitf1nch;mnch1]

alphaCh0 = mch0(1,2) / (mch0(1,2) + mch0(1,3))
sigmaCh0 = 1 / (mch0(1,2) + mch0(1,3))
biasCh0rl = mch0(1,1) / (mch0(1,2) + mch0(1,3))
biasCh0healthier = mch0(1,4) / (mch0(1,2) + mch0(1,3))

alphanCh0 = mnch0(1,2) / (mnch0(1,2) + mnch0(1,3))
sigmanCh0 = 1 / (mnch0(1,2) + mnch0(1,3))
biasnCh0rl = mnch0(1,1) / (mnch0(1,2) + mnch0(1,3))
biasnCh0healthier = mnch0(1,4) / (mnch0(1,2) + mnch0(1,3))

alphaCh1 = mch1(1,2) / (mch1(1,2) + mch1(1,3))
sigmaCh1 = 1 / (mch1(1,2) + mch1(1,3))
biasCh1rl = mch1(1,1) / (mch1(1,2) + mch1(1,3))
biasCh1healthier = mch1(1,4) / (mch1(1,2) + mch1(1,3))

alphanCh1 = mnch1(1,2) / (mnch1(1,2) + mnch1(1,3))
sigmanCh1 = 1 / (mnch1(1,2) + mnch1(1,3))
biasnCh1rl = mnch1(1,1) / (mnch1(1,2) + mnch1(1,3))
biasnCh1healthier = mnch1(1,4) / (mnch1(1,2) + mnch1(1,3))


para0found = [para0found;alphaCh0, biasCh0rl, sigmaCh0, biasCh0healthier]
para1found = [para1found;alphaCh1, biasCh1rl, sigmaCh1, biasCh1healthier]

para0found_n = [para0found_n; alphanCh0, biasnCh0rl, sigmanCh0, biasnCh0healthier]
para1found_n = [para1found_n; alphanCh1, biasnCh1rl, sigmanCh1, biasnCh1healthier]

hmch0 = glmfit([hrl0 rhealthier0], right0.*ch0, 'binomial','Link','probit')';
hmch1 = glmfit([hrl1 rhealthier1], right1.*ch1, 'binomial','Link','probit')';
hmnch0 = glmfit([hrl0 rhealthier0], right0.*nch0, 'binomial','Link','probit')';
hmnch1 = glmfit([hrl1 rhealthier1], right1.*nch1, 'binomial','Link','probit')';

hprobitf0ch = [hprobitf0ch;hmch0]
hprobitf1ch = [hprobitf1ch;hmch1]
hprobitf0nch = [hprobitf0nch;hmnch0]
hprobitf1nch = [hprobitf1nch;hmnch1]

tamch0 = glmfit([trl0 rhealthier0], right0.*ch0, 'binomial','Link','probit')';
tamch1 = glmfit([trl1 rhealthier1], right1.*ch1, 'binomial','Link','probit')';
tamnch0 = glmfit([trl0 rhealthier0], right0.*nch0, 'binomial','Link','probit')';
tamnch1 = glmfit([trl1 rhealthier1] , right1.*nch1, 'binomial','Link','probit')';

taprobitf0ch = [taprobitf0ch;tamch0]
taprobitf1ch = [taprobitf1ch;tamch1]
taprobitf0nch = [taprobitf0nch;tamnch0]
taprobitf1nch = [taprobitf1nch;tamnch1]

end


%explore probit regressions

%first full model with biasRL and bias H/T

mch=[]
mnch=[]

for o = 1:50
    
"-------SUBJECT---------two coeffs"  
o


"mch0"
mch0 =glmfit([hrl0all(:,o) trl0all(:,o) rhealthier0all(:,o)], right0all(:,o).*tch0(:,o), 'binomial','Link','probit')';
"mch1"
mch1 = glmfit([hrl1all(:,o) trl1all(:,o) rhealthier1all(:,o)], right1all(:,o).*tch1(:,o), 'binomial','Link','probit')';
"mnch0"
mnch0 = glmfit([hrl0all(:,o) trl0all(:,o) rhealthier0all(:,o)], right0all(:,o).*tnch0(:,o), 'binomial','Link','probit')';
"mnch1"
mnch1 = glmfit([hrl1all(:,o) trl1all(:,o) rhealthier1all(:,o)], right1all(:,o).*tnch1(:,o), 'binomial','Link','probit')';

"-------SUBJECT---------HD only"  
o


"mch0h"
mch0h =glmfit([hrl0all(:,o)  rhealthier0all(:,o)], right0all(:,o).*tch0(:,o), 'binomial','Link','probit')';
"mch1h"
mch1h = glmfit([hrl1all(:,o)  rhealthier1all(:,o)], right1all(:,o).*tch1(:,o), 'binomial','Link','probit')';
"mnch0h"
mnch0h = glmfit([hrl0all(:,o)  rhealthier0all(:,o)], right0all(:,o).*tnch0(:,o), 'binomial','Link','probit')';
"mnch1h"
mnch1h = glmfit([hrl1all(:,o)  rhealthier1all(:,o)], right1all(:,o).*tnch1(:,o), 'binomial','Link','probit')';

"-------SUBJECT---------TD only"  
o


"mch0t"
mch0t =glmfit([ trl0all(:,o) rhealthier0all(:,o)], right0all(:,o).*tch0(:,o), 'binomial','Link','probit')';
"mch1t"
mch1t = glmfit([ trl1all(:,o) rhealthier1all(:,o)], right1all(:,o).*tch1(:,o), 'binomial','Link','probit')';
"mnch0t"
mnch0t = glmfit([ trl0all(:,o) rhealthier0all(:,o)], right0all(:,o).*tnch0(:,o), 'binomial','Link','probit')';
"mnch1t"
mnch1t = glmfit([ trl1all(:,o) rhealthier1all(:,o)], right1all(:,o).*tnch1(:,o), 'binomial','Link','probit')';


%commands for checking the coeffs signs
%mch=[mch;mch0(1,2:3),mch1(1,2:3),mch0h(1,2),mch1h(1,2), mch0t(1,2),mch1t(1,2)]
%mnch=[mnch;mnch0(1,2:3),mnch1(1,2:3),mnch0h(1,2),mnch1h(1,2), mnch0t(1,2),mnch1t(1,2)]

end

%unidentifiable probit coeffs:
% 23 nch0, 21 nch0 nch1, 19 nch1, 18 nch0, 15 nch1, 14 nch1, 13 nch0, 12 nch1, 5 nch0,  
% 49 nch1, 48 nch0 nch1, 42 nch1, 40 nch0 nch1, 39 nch0, 37 nch1, 26 nch0 nch1, 
% 8 nch1: (full model unidentified with two large/unidentified positive coeffs in full model so can't remove one of them) 
%3 nch1 and 17 nch1: no non-negative coefficients (with at least one positive) found

%Use long dataset to run mixed probit (in Stata) to get group-level
%coefficients for the above subjects with unidentifiable probit coeffs for
%non-chall trials

%condition 0
 hrl_0 = .2800541  
 trl_0 = .183878   
 rhealthier_0 = .2042806  
 const_0 = -.0878851    

alphanCh0_gr = (hrl_0/trl_0) / (1 + hrl_0/trl_0)
sigmanCh0_gr = alphanCh0_gr / hrl_0
biasnCh0rl_gr = const_0 * (alphanCh0_gr / hrl_0)
biasnCh0healthier_gr = rhealthier_0 * (alphanCh0_gr / hrl_0)
 
%condition 1
 hrl_1 = .2071166
 trl_1 = .1516648  
 rhealthier_1 = .3265391  
 const_1 = .0050863 

 alphanCh1_gr = (hrl_1/trl_1) / (1 + hrl_1/trl_1)
sigmanCh1_gr = alphanCh1_gr / hrl_1
biasnCh1rl_gr = const_1 * (alphanCh1_gr / hrl_1)
biasnCh1healthier_gr = rhealthier_1 * (alphanCh1_gr / hrl_1)

for o = [3 5 8 12 13 14 15 17 18 19 21 23 26 37 39 40 42 48 49]
   
    para0found_n(o,1) = alphanCh0_gr
    para0found_n(o,3) = sigmanCh0_gr
    para0found_n(o,2) = biasnCh0rl_gr
    para0found_n(o,4) = biasnCh0healthier_gr
    
    para1found_n(o,1) = alphanCh1_gr
    para1found_n(o,3) = sigmanCh1_gr
    para1found_n(o,2) = biasnCh1rl_gr
    para1found_n(o,4) = biasnCh1healthier_gr 
  
end


%macc is now for cut-off for 0
macc = 0.0000000000001

%%SWITCH BETWEEN CHALL AND NCHALL where %FOR CHALL / %FOR NCHALL

id=[]

prob_oo = []
probk=[]

    allsigmaH0=[]
    allsigmaT0=[]
    allsigmaH1=[]
    allsigmaT1=[]
    allro0=[ ]
    allro1=[]
    
    allsigmaH0_n=[]
    allsigmaT0_n=[]
    allsigmaH1_n=[]
    allsigmaT1_n=[]
    allro0_n=[]
    allro1_n=[]
    
for o = [1:50]

    id=[id;o]
    
Phj = []
Ptj=[]
scprobj=[] 
spprobj=[]

    jsigmaH0=[]
    jsigmaT0=[]
    jsigmaH1=[]
    jsigmaT1=[]
    jro0=[]
    jro1=[]
    
    jsigmaH0_n=[]
    jsigmaT0_n=[]
    jsigmaH1_n=[]
    jsigmaT1_n=[]
    jro0_n=[]
    jro1_n=[]
    
    for j = 0:0.001:1
              
     
     %FOR CHALL
    %para0found_n(o,1) = j
    %para0found_n(o,3) = macc
    
    %para1found_n(o,1) = j
    %para1found_n(o,3) = macc
    
    %FOR NCHALL
    para0found(o,1) = j
    para0found(o,3) = macc
    
    para1found(o,1) = j
    para1found(o,3) = macc
    
%condition 0
if para0found(o,1) < macc & para0found_n(o,1) < macc
sigmaT0(o,1)=para0found(o,3)
sigmaT0_n(o,1)=para0found_n(o,3)
sigmaH0(o,1)= 10000
sigmaH0_n(o,1)= 10000
ro0(o,1)=0

elseif para0found(o,1) > 1 - macc & para0found_n(o,1) > 1 - macc
sigmaH0(o,1)=para0found(o,3)
sigmaH0_n(o,1)=para0found_n(o,3)
sigmaT0(o,1)= 10000
sigmaT0_n(o,1)= 10000
ro0(o,1)=0

elseif para0found(o,1) > 1 - macc & para0found_n(o,1) < macc
sigmaH0(o,1)=para0found(o,3)
sigmaH0_n(o,1)=10000
sigmaT0(o,1)= 10000
sigmaT0_n(o,1)= para0found_n(o,3)
ro0(o,1)=0

elseif para0found(o,1) < macc & para0found_n(o,1) > 1 - macc
sigmaH0(o,1)=10000
sigmaH0_n(o,1)=para0found_n(o,3)
sigmaT0(o,1)= para0found(o,3)
sigmaT0_n(o,1)= 10000
ro0(o,1)=0

else

%A means as in solution a) as in the text (Model section)

ro0A(o,1) = -sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1)))   
if ro0A(o,1)<0
ro0(o,1) = max(ro0A(o,1), 1/ro0A(o,1))
else
    ro0(o,1)=0
end
indA0(o,1) = 1 - abs( sign (ro0(o,1) - ro0A(o,1)) )  

sigmaH0f(o,1)=(para0found(o,3)/para0found(o,1))/sqrt(1-ro0(o,1)^2)
sigmaT0f(o,1)=(para0found(o,3)/(1-para0found(o,1)))/sqrt(1-ro0(o,1)^2)
sigmaH0f_n(o,1)=(para0found_n(o,3)/para0found_n(o,1))/sqrt(1-ro0(o,1)^2)
sigmaT0f_n(o,1)=(para0found_n(o,3)/(1-para0found_n(o,1)))/sqrt(1-ro0(o,1)^2)

if sigmaH0f(o,1)==inf
    sigmaH0f(o,1)=10000
else
sigmaH0f(o,1) = sigmaH0f(o,1)
end

if sigmaT0f(o,1)==inf
    sigmaT0f(o,1)=10000
else
sigmaT0f(o,1) = sigmaT0f(o,1)
end

if sigmaH0f_n(o,1)==inf
    sigmaH0f_n(o,1)=10000
else
sigmaH0f_n(o,1) = sigmaH0f_n(o,1)
end

if sigmaT0f_n(o,1)==inf
    sigmaT0f_n(o,1)=10000
else
sigmaT0f_n(o,1) = sigmaT0f_n(o,1)
end

if para0found(o,1) > 1 - macc
   fsigmaH0f(o,1) = macc
   fsigmaT0f(o,1) = 10000
elseif para0found(o,1) < macc
     fsigmaH0f(o,1) = 10000
     fsigmaT0f(o,1) = macc
else
fsigmaT0f(o,1)=sigmaH0f(o,1) * (abs(ro0(o,1)) / ((1-para0found(o,1))/para0found(o,1)))
fsigmaH0f(o,1)=sigmaT0f(o,1) * (abs(ro0(o,1)) / (para0found(o,1)/(1-para0found(o,1))))
end

if para0found_n(o,1) > 1 - macc
   fsigmaH0f_n(o,1) = macc
   fsigmaT0f_n(o,1) = 10000
elseif para0found_n(o,1) < macc
     fsigmaH0f_n(o,1) = 10000
     fsigmaT0f_n(o,1) = macc
else
fsigmaT0f_n(o,1)=sigmaH0f_n(o,1) * (abs(ro0(o,1)) / ((1-para0found_n(o,1))/para0found_n(o,1)))
fsigmaH0f_n(o,1)=sigmaT0f_n(o,1) * (abs(ro0(o,1)) / (para0found_n(o,1)/(1-para0found_n(o,1))))
end

sigmaH0(o,1)   = indA0(o,1)*fsigmaH0f(o,1)  + (1-indA0(o,1))*sigmaH0f(o,1)         %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaH0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaH0f(o,1)
sigmaH0_n(o,1) = indA0(o,1)*sigmaH0f_n(o,1) + (1-indA0(o,1))*fsigmaH0f_n(o,1)       %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaH0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaH0f_n(o,1)
sigmaT0(o,1)   = indA0(o,1)*sigmaT0f(o,1)   + (1-indA0(o,1))*fsigmaT0f(o,1)         % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaT0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaT0f(o,1)
sigmaT0_n(o,1) = indA0(o,1)*fsigmaT0f_n(o,1)  + (1-indA0(o,1))*sigmaT0f_n(o,1)          % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaT0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaT0f_n(o,1)

if sigmaH0(o,1)==0
     sigmaH0(o,1)=macc
else
sigmaH0(o,1)=sigmaH0(o,1)
end

if sigmaH0_n(o,1)==0
     sigmaH0_n(o,1)=macc
else
sigmaH0_n(o,1)=sigmaH0_n(o,1)
end

if sigmaT0(o,1)==0
     sigmaT0(o,1)=macc
else
sigmaT0(o,1)=sigmaT0(o,1)
end

if sigmaT0_n(o,1)==0
     sigmaT0_n(o,1)=macc
else
sigmaT0_n(o,1)=sigmaT0_n(o,1)
end



end

%condition 1
if para1found(o,1) < macc & para1found_n(o,1) < macc
sigmaT1(o,1)=para1found(o,3)
sigmaT1_n(o,1)=para1found_n(o,3)
sigmaH1(o,1)= 10000
sigmaH1_n(o,1)= 10000
ro1(o,1)=0

elseif para1found(o,1) > 1 - macc & para1found_n(o,1) > 1 - macc
sigmaH1(o,1)=para1found(o,3)
sigmaH1_n(o,1)=para1found_n(o,3)
sigmaT1(o,1)= 10000
sigmaT1_n(o,1)= 10000
ro1(o,1)=0

elseif para1found(o,1) > 1 - macc & para1found_n(o,1) < macc
sigmaH1(o,1)=para1found(o,3)
sigmaH1_n(o,1)=10000
sigmaT1(o,1)= 10000
sigmaT1_n(o,1)= para1found_n(o,3)
ro1(o,1)=0

elseif para1found(o,1) < macc & para1found_n(o,1) > 1 - macc
sigmaH1(o,1)=10000
sigmaH1_n(o,1)=para1found_n(o,3)
sigmaT1(o,1)= para1found(o,3)
sigmaT1_n(o,1)= 10000
ro1(o,1)=0

else
    
%ro1A(o,1) = -sqrt( (para1found_n(o,1)/(1-para1found_n(o,1))) * ((1-para1found(o,1))/para1found(o,1)))   
ro1A(o,1) = -sqrt( (para1found(o,1)/(1-para1found(o,1))) * ((1-para1found_n(o,1))/para1found_n(o,1)))   
if ro1A(o,1)<0
ro1(o,1) = max(ro1A(o,1), 1/ro1A(o,1))
else
    ro1(o,1)=0
end
indA1(o,1) = 1 - abs( sign (ro1(o,1) - ro1A(o,1)) ) 

sigmaH1f(o,1)=(para1found(o,3)/para1found(o,1))/sqrt(1-ro1(o,1)^2)
sigmaT1f(o,1)=(para1found(o,3)/(1-para1found(o,1)))/sqrt(1-ro1(o,1)^2)
sigmaH1f_n(o,1)=(para1found_n(o,3)/para1found_n(o,1))/sqrt(1-ro1(o,1)^2)
sigmaT1f_n(o,1)=(para1found_n(o,3)/(1-para1found_n(o,1)))/sqrt(1-ro1(o,1)^2)

if sigmaH1f(o,1)==inf
    sigmaH1f(o,1)=10000
else
sigmaH1f(o,1) = sigmaH1f(o,1)
end

if sigmaT1f(o,1)==inf
    sigmaT1f(o,1)=10000
else
sigmaT1f(o,1) = sigmaT1f(o,1)
end

if sigmaH1f_n(o,1)==inf
    sigmaH1f_n(o,1)=10000
else
sigmaH1f_n(o,1) = sigmaH1f_n(o,1)
end

if sigmaT1f_n(o,1)==inf
    sigmaT1f_n(o,1)=10000
else
sigmaT1f_n(o,1) = sigmaT1f_n(o,1)
end

if para1found(o,1) > 1 - macc
   fsigmaH1f(o,1) = macc
   fsigmaT1f(o,1) = 10000
elseif para1found(o,1) < macc
     fsigmaH1f(o,1) = 10000
     fsigmaT1f(o,1) = macc
else
fsigmaT1f(o,1)=sigmaH1f(o,1) * (abs(ro1(o,1)) / ((1-para1found(o,1))/para1found(o,1)))
fsigmaH1f(o,1)=sigmaT1f(o,1) * (abs(ro1(o,1)) / (para1found(o,1)/(1-para1found(o,1))))
end

if para1found_n(o,1) > 1 - macc
   fsigmaH1f_n(o,1) = macc
   fsigmaT1f_n(o,1) = 10000
elseif para1found_n(o,1) < macc
     fsigmaH1f_n(o,1) = 10000
     fsigmaT1f_n(o,1) = macc
else
fsigmaT1f_n(o,1)=sigmaH1f_n(o,1) * (abs(ro1(o,1))/ ((1-para1found_n(o,1))/para1found_n(o,1)))
fsigmaH1f_n(o,1)=sigmaT1f_n(o,1) * (abs(ro1(o,1)) / (para1found_n(o,1)/(1-para1found_n(o,1))))
end


sigmaH1(o,1)   = indA1(o,1)*fsigmaH1f(o,1)  + (1-indA1(o,1))*sigmaH1f(o,1)         %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaH0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaH0f(o,1)
sigmaH1_n(o,1) = indA1(o,1)*sigmaH1f_n(o,1) + (1-indA1(o,1))*fsigmaH1f_n(o,1)       %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaH0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaH0f_n(o,1)
sigmaT1(o,1)   = indA1(o,1)*sigmaT1f(o,1)   + (1-indA1(o,1))*fsigmaT1f(o,1)         % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaT0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaT0f(o,1)
sigmaT1_n(o,1) = indA1(o,1)*fsigmaT1f_n(o,1)  + (1-indA1(o,1))*sigmaT1f_n(o,1)          % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaT0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaT0f_n(o,1)

if sigmaH1(o,1)==0
     sigmaH1(o,1)=macc
else
sigmaH1(o,1)=sigmaH1(o,1)
end

if sigmaH1_n(o,1)==0
     sigmaH1_n(o,1)=macc
else
sigmaH1_n(o,1)=sigmaH1_n(o,1)
end

if sigmaT1(o,1)==0
     sigmaT1(o,1)=macc
else
sigmaT1(o,1)=sigmaT1(o,1)
end

if sigmaT1_n(o,1)==0
     sigmaT1_n(o,1)=macc
else
sigmaT1_n(o,1)=sigmaT1_n(o,1)
end



end





 %Parameters_challenging = ["alpha0", "ro0", "sigmaH0", "sigmaT0", "sigma0","bias0rl", "bias0healthier", "alpha1","ro1", "sigmaH1", "sigmaT1", "sigma1", "bias1rl", "bias1healthier"; para0found(:,1), ro0, sigmaH0, sigmaT0, para0found(:,3),para0found(:,2), para0found(:,4),para1found(:,1), ro1, sigmaH1, sigmaT1, para1found(:,3),para1found(:,2), para1found(:,4)]
 %Parameters_nonchallenging = ["alpha0n", "ro0n", "sigmaH0n", "sigmaT0n", "sigma0n","bias0nrl", "bias0nhealthier", "alpha1n","ro1n", "sigmaH1n", "sigmaT1n", "sigma1n", "bias1nrl", "bias1nhealthier"; para0found_n(:,1), ro0, sigmaH0_n, sigmaT0_n, para0found_n(:,3),para0found_n(:,2), para0found_n(:,4),para1found_n(:,1), ro1, sigmaH1_n, sigmaT1_n, para1found_n(:,3),para1found_n(:,2), para1found_n(:,4)]



       
    Ph0_o = []  
    Pt0_o = [] 
    Ph1_o = []  
    Pt1_o = []
    spprob0_o = []
    scprob0_o = []
    spprob1_o = []
    scprob1_o = []
         
        
     Ph0_ou = []  
    Pt0_ou = [] 
    Ph1_ou = []  
    Pt1_ou = []
    spprob0_ou = []
    scprob0_ou = []
    spprob1_ou = []
    scprob1_ou = []
    
     for i = [1:90]
       
         % differences for chosen minus unchosen item
        mu0 = [right0all(i,o)*hrl0all(i,o) + (1-right0all(i,o))*(-hrl0all(i,o)),right0all(i,o)*trl0all(i,o) + (1-right0all(i,o))*(-trl0all(i,o))]
        mu1 = [right1all(i,o)*hrl1all(i,o) + (1-right1all(i,o))*(-hrl1all(i,o)),right1all(i,o)*trl1all(i,o) + (1-right1all(i,o))*(-trl1all(i,o))]
                       
        %FOR CHALL
        %sigma0 =  [sigmaH0(o,1), ro0(o,1); ro0(o,1),sigmaT0(o,1)]
        %sigma1 =  [sigmaH1(o,1), ro1(o,1); ro1(o,1),sigmaT1(o,1)]
        
        %FOR NCHALL AND add o~=47 for sigma "if.." [sigma impossible for
        %at least some alphaQ intervals
        sigma0 =  [sigmaH0_n(o,1) , ro0(o,1); ro0(o,1) ,sigmaT0_n(o,1)]
        sigma1 =  [sigmaH1_n(o,1) , ro1(o,1); ro1(o,1) ,sigmaT1_n(o,1)]
        
        if sigma0(1,1)>0&sigma0(1,2)>-1.01&sigma0(2,1)>-1.01&sigma0(2,2)>0&sigma0(1,1)<inf&sigma0(2,2)<inf&o~=47
        Ph0_o = [Ph0_o; mvncdf([0 -inf],[inf inf],mu0,sigma0)]
        Pt0_o = [Pt0_o; mvncdf([-inf 0],[inf inf],mu0,sigma0)]
        spprob0_o = [spprob0_o; mvncdf([0 0],[inf inf],mu0,sigma0) / (mvncdf([0 0],[inf inf],mu0,sigma0) + mvncdf([-inf -inf],[0 0],mu0,sigma0))]
        scprob0_o = [scprob0_o; mvncdf([0 -inf],[inf 0],mu0,sigma0) / (mvncdf([0 -inf],[inf 0],mu0,sigma0) + mvncdf([-inf 0],[0 inf],mu0,sigma0))]
                  
                  
                                 
        else
            Ph0_o = [Ph0_o; NaN]
            Pt0_o = [Pt0_o; NaN]
            spprob0_o = [spprob0_o; NaN]
            scprob0_o = [scprob0_o; NaN]
          
        end
        
        
        if sigma1(1,1)>0&sigma1(1,2)>-1.01&sigma1(2,1)>-1.01&sigma1(2,2)>0&sigma1(1,1)<inf&sigma1(2,2)<inf&o~=47
         Ph1_o = [Ph1_o; mvncdf([0 -inf],[inf inf],mu1,sigma1)]
        Pt1_o =  [Pt1_o; mvncdf([-inf 0],[inf inf],mu1,sigma1)]
        spprob1_o = [spprob1_o; mvncdf([0 0],[inf inf],mu1,sigma1) / (mvncdf([0 0],[inf inf],mu1,sigma1) + mvncdf([-inf -inf],[0 0],mu1,sigma1))]
        scprob1_o = [scprob1_o; mvncdf([0 -inf],[inf 0],mu1,sigma1) / (mvncdf([0 -inf],[inf 0],mu1,sigma1) + mvncdf([-inf 0],[0 inf],mu1,sigma1))]
               
       
                               
            
        else
            Ph1_o = [Ph1_o; NaN]
            Pt1_o = [Pt1_o; NaN]
            spprob1_o = [spprob1_o; NaN]
            scprob1_o = [scprob1_o; NaN]
            
        end
        
       
   
        
        
     end
     
     
   
     
          Ph = [Ph0_o; Ph1_o]
          Pt = [Pt0_o; Pt1_o]
          spprob = [spprob0_o; spprob1_o]
          scprob = [scprob0_o; scprob1_o]
          
          
   chall = [tchall0(:,o);tchall1(:,o)]
   
   for i=1:180
              if chall(i,1)==0 %FOR CHALL PUT 1, FOR NCHALL PUT 0
                  Ph(i,1)=Ph(i,1)
                  Pt(i,1)=Pt(i,1)
                  spprob(i,1)=spprob(i,1)
                  scprob(i,1)=scprob(i,1)
              else
                  Ph(i,1)=NaN
                  Pt(i,1)=NaN
                  spprob(i,1)=NaN
                  scprob(i,1)=NaN
              end
          end
          
   
   
   
   cond = [zeros(90,1); ones(90,1)]
   prob_o = [o.*ones(180,1), cond, chall] 
   
Phj = [Phj, Ph]
Ptj = [Ptj, Pt]
scprobj = [scprobj, scprob]
spprobj = [spprobj, spprob]


    jsigmaH0=[jsigmaH0, sigmaH0(o,1)]
    jsigmaT0=[jsigmaT0, sigmaT0(o,1)]
    jsigmaH1=[jsigmaH1, sigmaH1(o,1)]
    jsigmaT1=[jsigmaT1, sigmaT1(o,1)]
    jro0=[jro0, ro0(o,1)]
    jro1=[jro1, ro1(o,1)]

    
    jsigmaH0_n=[jsigmaH0_n, sigmaH0_n(o,1)]
    jsigmaT0_n=[jsigmaT0_n, sigmaT0_n(o,1)]
    jsigmaH1_n=[jsigmaH1_n, sigmaH1_n(o,1)]
    jsigmaT1_n=[jsigmaT1_n, sigmaT1_n(o,1)]
    jro0_n=[jro0_n,ro0(o,1) ]
    jro1_n=[jro1_n, ro1(o,1)]

    end  
    probk=[probk;nanmean(Phj,2), nanmean(Ptj,2),nanmean(scprobj,2),nanmean(spprobj,2)]
    prob_oo = [prob_oo; prob_o]
    
    
    allsigmaH0=[allsigmaH0; jsigmaH0]
    allsigmaT0=[allsigmaT0; jsigmaT0]
    allsigmaH1=[allsigmaH1; jsigmaH1]
    allsigmaT1=[allsigmaT1; jsigmaT1]
    allro0=[allro0;jro0 ]
    allro1=[allro1; jro1]
    
    
    allsigmaH0_n=[allsigmaH0_n; jsigmaH0_n]
    allsigmaT0_n=[allsigmaT0_n; jsigmaT0_n]
    allsigmaH1_n=[allsigmaH1_n; jsigmaH1_n]
    allsigmaT1_n=[allsigmaT1_n; jsigmaT1_n]
    allro0_n=[allro0_n;jro0_n ]
    allro1_n=[allro1_n; jro1_n]   
      
end

tprob = ["subject", "condition", "chall", "Ph","Pt", "scprob", "spprob"; prob_oo, probk]  
 
 
   msigmaH0=mean(allsigmaH0,2)
    msigmaT0=mean(allsigmaT0,2)
    msigmaH1=mean(allsigmaH1,2)
    msigmaT1=mean(allsigmaT1,2)
    mro0=mean(allro0,2)
    mro1=mean(allro1,2)

        
msigmaH0_n=mean(allsigmaH0_n,2)
    msigmaT0_n=mean(allsigmaT0_n,2)
    msigmaH1_n=mean(allsigmaH1_n,2)
    msigmaT1_n=mean(allsigmaT1_n,2)
    mro0_n=mean(allro0_n,2)
    mro1_n=mean(allro1_n,2)
    
   
Parameters_challenging = ["alpha0", "ro0", "sigmaH0", "sigmaT0", "sigma0","bias0rl", "bias0healthier", "alpha1","ro1", "sigmaH1", "sigmaT1", "sigma1", "bias1rl", "bias1healthier"; para0found(:,1), mro0, msigmaH0, msigmaT0, para0found(:,3),para0found(:,2), para0found(:,4),para1found(:,1), mro1, msigmaH1, msigmaT1, para1found(:,3),para1found(:,2), para1found(:,4)]

Parameters_nonchallenging = ["alpha0n", "ro0n", "sigmaH0n", "sigmaT0n", "sigma0n","bias0nrl", "bias0nhealthier", "alpha1n","ro1n", "sigmaH1n", "sigmaT1n", "sigma1n", "bias1nrl", "bias1nhealthier"; para0found_n(:,1), mro0_n, msigmaH0_n, msigmaT0_n, para0found_n(:,3),para0found_n(:,2), para0found_n(:,4),para1found_n(:,1), mro1_n, msigmaH1_n, msigmaT1_n, para1found_n(:,3),para1found_n(:,2), para1found_n(:,4)]


save('Approach1_NCHALL_1000')